{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stop Hotspot Detection\n",
    "\n",
    "<img align=\"right\" src=\"https://anitagraser.github.io/movingpandas/assets/img/movingpandas.png\">\n",
    "\n",
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/anitagraser/movingpandas-examples/main?filepath=3-analysis-examples/stop-hotspots.ipynb)\n",
    "\n",
    "This demo shows how to detect significant hotspots of urban activity. It combines MovingPandas, OSMnx, and Spaghetti and was inspired by:\n",
    "\n",
    "* [Network-constrained spatial autocorrelation by James D. Gaboardi](https://github.com/pysal/spaghetti/blob/8dfec40812a9824e575822a1587b47f26fc45313/notebooks/network-spatial-autocorrelation.ipynb)\n",
    "* [Simplify network topology and consolidate intersections by Geoff Boeing](https://github.com/gboeing/osmnx-examples/blob/64e104f8e3e719c23c640172c2f18ba7b46a020d/notebooks/04-simplify-graph-consolidate-nodes.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "from osmnx import utils_graph\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import movingpandas as mpd\n",
    "from datetime import timedelta\n",
    "import spaghetti\n",
    "import esda\n",
    "from splot.esda import moran_scatterplot, lisa_cluster, plot_moran, plot_local_autocorrelation\n",
    "\n",
    "ox.config(use_cache=True, log_console=True)\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding stops in movement trajectories\n",
    "\n",
    "This also serves as a stress test for MovingPandas. The sample dataset contains 190k points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks_gdf = gpd.read_file('../data/geolife_sample.gpkg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks_gdf['t'] = pd.to_datetime(tracks_gdf['t'])\n",
    "tracks_gdf = tracks_gdf.set_index('t').tz_localize(None)\n",
    "traj_collection = mpd.TrajectoryCollection(tracks_gdf, 'trajectory')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "stop_pts = mpd.TrajectoryStopDetector(traj_collection)\\\n",
    "   .get_stop_points(min_duration=timedelta(seconds=180), max_diameter=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(stop_pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minx, miny, maxx, maxy = stop_pts.geometry.total_bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing the network data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "buffer = 0.005\n",
    "G = ox.graph_from_bbox(maxy+buffer, miny-buffer, minx-buffer, maxx+buffer, network_type='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_proj = ox.project_graph(G)  # projects to suitable UTM CRS\n",
    "G2 = ox.consolidate_intersections(G_proj, rebuild_graph=True, tolerance=15, dead_ends=False)\n",
    "G2 = ox.project_graph(G2, to_crs='EPSG:4326')\n",
    "len(G2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph_gdf = gpd.GeoDataFrame(utils_graph.graph_to_gdfs(G2, nodes=False)[\"geometry\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph_gdf.plot(figsize=(9,9), color='grey')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's worth noting that there's an area in the east (in the Temple of Heaven park) with a lot of stops but no network coverage. This happend because these areas are modelled as open spaces (highway=pedestrian, area=true) in OSM (https://www.openstreetmap.org/way/29182811)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = graph_gdf.plot(figsize=(9,9), color='grey', zorder=0)\n",
    "stop_pts.plot(ax=ax, color='red', zorder=3, alpha=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assigning stops to network edges "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ntw = spaghetti.Network(in_data=graph_gdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_name = \"stops\"  # point pattern name\n",
    "ntw.snapobservations(stop_pts, pp_name, attribute=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating local statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_moran(net, pp_name, w):\n",
    "    \"\"\"Calculate a Moran's I statistic based on network arcs.\"\"\"\n",
    "    # Compute the counts\n",
    "    pointpat = net.pointpatterns[pp_name]\n",
    "    counts = net.count_per_link(pointpat.obs_to_arc, graph=False)\n",
    "    # Build the y vector\n",
    "    arcs = w.neighbors.keys()\n",
    "    y = [counts[a] if a in counts.keys() else 0. for i, a in enumerate(arcs)]\n",
    "    # Moran's I\n",
    "    moran = esda.moran.Moran(y, w, permutations=99)\n",
    "    return moran, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moran_ntwwn, yaxis_ntwwn = calc_moran(ntw, pp_name, ntw.w_network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, arc_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)\n",
    "arc_df['n'] = yaxis_ntwwn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = arc_df.plot(column='n', figsize=(12,9), cmap='Reds', legend=True, zorder=3, lw=3, vmax=10)\n",
    "stop_pts.plot(ax=ax, color='black', alpha=0.2, zorder=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moran_loc_ntwwn = esda.moran.Moran_Local(yaxis_ntwwn, ntw.w_network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = lisa_cluster(moran_loc_ntwwn, arc_df, p=0.05, figsize=(9,9), lw=3, zorder=3)\n",
    "stop_pts.plot(ax=ax, zorder=1, alpha=.2, color=\"black\", markersize=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continue exploring MovingPandas\n",
    "\n",
    "\n",
    "1. [Bird migration analysis](1-bird-migration.ipynb)\n",
    "1. [Ship data analysis](2-ship-data.ipynb)\n",
    "1. [Horse collar data exploration](3-horse-collar.ipynb)\n",
    "1. [Stop hotspot detection](4-stop-hotspots.ipynb)\n",
    "1. [OSM traces](5-osm-traces.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
